Holstein polaron: the effect of multiple phonon modes 
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We generalize the Momentum Average approximations MA^ and MA^ to study the effects of 
coupling to multiple optical phonons on the properties of a Holstein polaron. As for a single phonon 
mode, these approximations are numerically very efficient. They become exact for very weak or very 
strong couplings, and are highly accurate in the intermediate regimes, e.g. the spectral weights obey 
exactly the first six, respectively eight, sum rules. Our results show that the effect on ground-state 
properties is cumulative in nature. In particular, if the effective coupling to one mode is much larger 
than to the others, this mode effectively determines the GS properties. However, even very weak 
coupling to a second phonon mode has important non-perturbational effects on the higher energy 
spectrum, in particular on the dispersion and the phonon statistics of the polaron band. 

PACS numbers: 71.38.-k, 72.10.Di, 63.20.Kr 



The coupling of electrons to phonons is a widely stud- 
ied problem, because it leads to many interesting phe- 
nomena such as conventional superconductivity or the 
formation of polarons (composite objects comprised of 
an electron and the surrounding phonon cloud), impor- 
tant in several classes of materials. As a recent example, 
results from angle-resolved photoemission spectroscopy^ 
(ARPES) have lead to new discussions about possible 
polaronic effects in high-temperature superconductors P 

Most theoretical studies of polaron properties are of 
the Holstein model with a single optical phonon modep^ 
even though complex materials have many optical and 
acoustic phonons. The reason is that usually there is one 
optical mode to which the coupling is strongest, and one 
assumes that the effects of the other modes are perturba- 
tionally small. Also, the efficiency of various numerical 
methods^ such as exact diagonalization, diagrammatic 
Monte Carlo and variational methods employed for ob- 
taining results in the intermediate coupling regime, where 
no exact solutions are known, suffers when the Hilbert 
space is enlarged by addition of multiple phonon modes. 

Recently, the so-called Momentum Average (MA) an- 
alytical approximation^ has been shown to be highly 
accurate over most of the parameter space of the Hol- 
stein polaron problem, while requiring a numerically triv- 
ial effort irrespective of the dimensionality of the problem 
or the strength of the coupling. Moreover, its accuracy 
can be systematically improved. 7 Such fast but accurate 
methods are useful for a quick survey of polaron prop- 
erties in various regimes, which can then be followed by 
quantitatively more accurate, but significantly more time 
and resource consuming numerical simulations. 

In this Rapid Communication we show how these ap- 
proximations can be generalized to deal with multiple 
phonon modes, without loss of accuracy when compared 
to the single-mode results. This allows us to study easily 
the effects of multiple phonon modes on the polaron prop- 
erties. As we show below, while for ground-state prop- 
erties these effects are rather trivial, the higher-energy 
spectrum is significantly modified by additional phonon 
modes, even if coupling to them is perturbationally small. 



The generalized Holstein Hamiltonian 8 of interest is: 

a,k,q 

The first term describes a free electron on a d-dimensional 
lattice with N sites (its spin is irrelevant, thus the spin in- 
dex is dropped), the second describes the optical phonon 
modes, and the third describes the electron-phonons cou- 
plings. Momenta sums are over the Brillouin zone. 
We focus on finding the polaron's Green's function!^ 

G(k, W ) = (0|c k GH4|0), (2) 

where G(uo) = [uo — H+ir]] -1 is the usual resolvent (% = 1) 
and |0) is the vacuum. The poles of this Green's function 
mark the polaron spectrum, and ground-state (GS) ener- 
gies, effective masses, quasiparticle (qp) weights and aver- 
age phonon numbers can then be calculated as discussed 
in Ref. 6. The spectral weight A(k,(j) = — ^G(k, uo) can 
also be directly compared against ARPES results. 

To simplify notation, we first assume that there are 
only two phonon modes, and rename their operators as 
6 q (mode 1) and Bq (mode 2). The generalization to 
more phonon modes is discussed below. We use repeat- 
edly Dyson's identity G(uo) = Go(uo) + G(uo)VGo(uo), 
where Go(uo) = [uo — Ho + if]] -1 corresponds to the non- 
interacting Hamiltonian, to generate an infinite system of 
coupled equations involving G(k, uo) and the generalized 
Green's functions F nm (k, qi, . . . , q n ; Qi, . . . , Q m ; uo) = 
(0\c^G(uj)cl T bl • • • ...B^ \0). Here k T = k - 

q T -Qr and q T = Y%=i Qt = Qj Arguments 

identical to those of Ref. [7] show that all functions F nm 
are proportional to G(k, uo). It is thus more convenient to 
work with the rescaled functions / nm (k, {q}, {Q}, uo) — 
A/^ n+m )/ 2 F nm (k, {q}, {Q},cj)/G(k,cj), where we use the 
shorthand notation {q} = q 1? . . . , q n , etc. 

The solution has the standard form 

G(k,<j) = ^-e k -S(k,cj) + ^] _1 (3) 
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where the exact self-energy is given by 

EM = |E /io (k, qi ,u>) + 1 J2 fa ( k > Qi > w ) W 
qi Qi 

I 



/nm({q},{Q}) = 



where the dependence on k, uj is implicitly assumed for 
all fnm-, Go(k,c<;) = (uj — ek + ^) _1 is the free propagator, 
and /oo = 1 by definition. These are the generalization 
of the equivalent single- mode equations of Ref. 7. 

As for the single-mode problem, the MA^ ^ approxi- 
mation is obtained by replacing in the r.h.s of Eqs. ([5| 

Go(kx, ou — nQi — ra^) ~^ 9o(^ — — 1T1Q2) = 9nm(^) 
where the momenta averages 

k 

are simple known functions.^ In terms of the functions 

^nm(^) = /nm(k,{q},{Q},^), (6) 

{q},{Q} 

the momentum-independent MA^ ^ self-energy is: 

^MAW = giFlo{v) + fl^OlM, 

while the recurrence relations ([5| take the simpler form 
(^)[n#iJF n _ 1?m (o;) + mg 2 J r n ,m-i(^) + 

tfl^n+l.mM +^n,m+lH] (of COUrse, J^ = !)• 

Such recursive equations were previously solved in a 
different context by Cini et a/., 10 but their solution can- 
not be generalized to more than two phonon modes, nor 
to MA^ 1 ) or higher levels (see below). We have found an 
alternative solution without these shortcomings. First, 
we rewrite these recurrence relations in matrix form: 

V k = A k V k -i+B k Vk+i, (7) 

where the vector V k = (J^.o; Fk-1,1] • • • ; Fi,k-i', ^o,k) T 
contains all k + 1 functions corresponding to a total 
of k phonons. is a matrix of size k + 1 x k with 
the only non-zero elements (A k )i,i = (k — i)gi(jk-i,i(w) 
and (Ak)i+i,i = (i + Vi = 0,fe - 1. 

Similarly, 1^ is a matrix of size fc + 1 x k + 2 with 
the only non-zero elements (B k ) i}i = gig k - iyi (uj) and 
= 929k-i 1 i(u), = 0,fc. Dependence on uj 
is again implicitly assumed everywhere. 



In terms of the new sets {q}^ = qi , . . . , 1 , q^ + i , . . . , q n 
and {q} n +i = qi, . . . , q n , q n +i, the functions f nm are 
the solutions of the following recurrence relations: 



(5) 
I 

The solution is V k = M k V k -i with Vq = (1), where 

M k = l - A k . (8) 

1 — B k A k+1 

1 — B k +i- A k +2 

is a continued fraction of matrices of increasing size. The 
self energy is S MA ( 0) (uj) = (gi,g 2 )Vi = (g 1 ,g 2 )M 1 . The 
continued fractions become convergent if truncated at 
levels N « gf/^i + g 2 /Q 2 , i- e - wnen one keeps contribu- 
tions from f nm corresponding to expected average num- 
bers of phonons in the cloud. The generalization to more 
phonon modes is straightforward. V k again contains all 
Green's functions with a fixed number of phonons, and 
the interaction links it only to V k ±\. The matrices A k 
and B k have non-vanishing elements only on a number 
of diagonals equal to the number of modes. The dimen- 
sion of these matrices increases now faster with increasing 
fc, but it is still much less severe than the corresponding 
increase in numerical simulations. Also, note that the 
MA calculation is equally simple in any dimension, the 
only change being in the expression used for go(u)W 

The analysis of the diagrammatic and variational 
meaning of MA^ ^ and the sum rule s it obeys, is iden- 
tical to that for the single-mode casej^^l anc l we d not 
repeat it. It again proves its accuracy over the entire 
parameter space, as long as Vti/t > 0.1, Vi. However, 
MA(°) fail s to c orrectly predict the polaron+one phonon 
continuum.^^ To remedy this, we use or a higher 

level approximation. 7 In MA^\ Eqs. ^ for / i and fio 
are left unchanged, and the momentum average is made 
only for f nrn with n + m > 2. As shown in Ref. [7J MA^ 1 ) 
correctly predicts the polaron+one phonon continuum, 
besides giving small improvements in the accuracy of var- 
ious other quantities, e.g. the number of exactly satisfied 
sum rules increases from six to eight. The derivation of 
S MA (i) (uj) follows identically that in Ref. 7, with the 
only difference that continued fractions there correspond 
to continued fractions of matrices [like in Eq. ^] here. 
We will present the details of this straightforward deriva- 
tion elsewhere, 12 instead focusing here on results. 



Go(kT,a; — nQi — mfl 2 ) 



91 fn-l,m({q}i, {Q}) + #2 Y fn 
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Y /n+l,m({q}n+l,{Q}) + Y /n,m+l ({q} , {Q}m+1 ) 
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All results shown are for two phonon modes and d = 1, 
which suffices to uncover the essential new physics. Re- 
sults for higher d and more phonon modes will be pre- 
sented elsewhere.^ We begin by discussing GS prop- 
erties such as the energy Eqs an d effective mass ra*, 
shown in Fig. [I] as functions of the effective couplings 
Xi = g? /(2dt£li), i = 1,2. Note that the qp weight is 
Zo = m/m*, where ra is the bare electron massP The 
"equipotential" lines drawn show that Eqs is we U de- 
scribed as a function of only A e fr = A^, whereas ra* 
is a function of ^2 i \i/O i . If Oi = O for all modes, 
these can be proved to be exact results.^ In the strong 
coupling limit one also expects Eqs — ~^2i9i A^i? 
InZo ex — Xi/fli, supporting the same conclusion. To 
a good extent the only effect of having ^ 2 is to 
change the slope of the ra* "equipotentials" from the 45° 
found when = although some slight deviations 
from linearity are also seen in Eqs when either A^ <C 1. 

We conclude that GS properties can be well under- 
stood in cumulative terms, for instance the energy is that 
of a polaron coupled to a single phonon with A e ff. The 
crossover from large to small-polaron behavior is there- 
fore expected when A e ff « lP=l As a result, it is possible 
to have small-polaron behavior even if each individual 
phonon mode is weakly coupled to the electron (A^ < 1). 
However, in cases where one mode (say, mode 1) is indeed 
much more strongly coupled than all others, Ai ^> A^, 
i = 2, . . . , then \ e ff ~ Ai and one can, to a good extent, 
ignore the small cumulative effect from the other modes. 

This conclusion, however, does not generally hold for 
higher-energy properties, as we show now. In Fig. [2] we 
plot the d = 1 spectral function A(k 1 uo) vs. k and uj. 
The effective coupling to the first mode, of frequency 
fli/t = 0.7, is kept constant to a fairly low value Ai = 0.4. 
In panel (a), coupling to the second mode, of energy 
2 /t = 0.3, is zero, so this is effectively a one-phonon 
problem. As expected, at low energies we see the po- 
laron band, of width 0\, followed above Eqs + ^i by 
the polaron+one-phonon continuum, and other features 



FIG. 2: (color online) A(k,w) from MA (1) for d = 1, and for 
0.7t, 2 0.3t, Ai 0.4 and a) A 2 0.0, b) A 2 0.2, 
c) A 2 = 0.6. 



at even higher energies ! 3 * 4 * 7 * 

Addition of even very weak coupling to a second 
phonon mode changes things considerably, as shown in 
panel (b) for A2 = 0.2. If 2 < Oi (as chosen here), the 
polaron band width is changed to ^2, even though the 
GS energy and effective mass are not much affected (see 
previous discussion). This significant change is not so 
surprising if one considers the origin of the polaron+one 
phonon continuum: it corresponds to states where one 
phonon is created far from the polaron. As a result, they 
interact little and the total energy is just the sum of 
the two. If there are several phonon modes, the contin- 
uum will be defined by the mode with the lowest fre- 
quency l^min, irrespective of whether this is the mode 
most strongly coupled to the electron or not. Because of 
this, the polaron band cannot be wider than ^ m i n . 

This interpretation is confirmed by phonon statis- 
tics, shown in Fig. [3j Here we plot average numbers 
Ni and N 2 of phonons of either type in the polaron 
cloud, as a function of the polaron momentum k. T hese 
are calculated using the Hellman-Feynman theoremP^' 
Again, the coupling to the first mode is kept constant 
at Ai = 0.4. If A2 = ( + symbols), we see that N\ 
increases from a small value at k — to just above 1 for 
k > 7r/2. This shows, as expected, that while around 
k = the large-polaron is essentially similar to a free 
electron, for k > tt/2 the largest contribution to the po- 
laron comes from electron+one phonon states. Of course, 
A^2 = in this case. As A2 is turned on but is still 
small (A2 = 0.1,0.3, x symbols) there is little change 
near k — 0, however the changes at higher momenta are 
dramatic: N\ decreases by 1 whereas N 2 increases by 1 
(note the different scales). This confirms that it is now 
the second type of phonon that controls the nature of 
the polaron at large k values, even though A2 < Ai. 
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FIG. 3: (color online) Number of phonons of either type in the 
polaron cloud vs. A2, for Qi = 0.7£, ^2 = 0.3£ and Ai = 0.4. 
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FIG. 4: (color online) A(k = 0, 00) vs. 00 and A 2 , for f2i = 0.7t, 
Q 2 = O.St and Ai = OAt. 



Once A2 > Ai, the second phonon completely domi- 
nates behavior. For A2 = 0.6 (square symbols) we have 
A e ff = Ai + A2 = 1, and the polaron is in the crossover 
regime, whereas for A2 = 1.0,A e ff = 1.4 the polaron is 
firmly in the small-polaron regime. In the latter case, 
we expect to see less and less k dependence, since the 
polaron cloud becomes limited to the site on which the 



electron resides. This behavior is indeed observed for 7V 2 , 
but also for N\. This is because once the polaron is lo- 
calized at a site (because of strong coupling to mode 2) it 
will automatically also shift the equilibrium position for 
mode 1 phonons at that site, resulting in the creation of 
a finite number of such bare phonons. 

If there is only one phonon mode coupled to the elec- 
tron, one can use ARPES data to extract its frequency 
Q from the location of the "discontinuity" in the disper- 
sion, at weak coupling, or from average distance between 
eigenstates at strong couplings, where the Lang-Firsov 
spectrum appears.^ The effective mass determines A, so 
g can also be extracted. Our results show that this pro- 
cedure is usually wrong in the case of multiple phonon 
modes, even if one expects coupling to one of them to 
be dominant. It only works if this particular mode also 
happens to have the lowest frequency, else one will un- 
derestimate its £1 and overestimate g. This point is 
clearly demonstrated in Fig.[4j which shows that the ^ m in 
phonon defines the location of the low-energy (not GS) 
features, irrespective of its coupling. Of course, one may 
hope to see the continua due to the other phonon modes 
(see also Fig. [5]c) and thus be able to identify their fre- 
quencies. This is probably unlikely, due to broadening in 
real data (note that temperature dependence would also 
be determined by the ^ m i n phonon, not by the dominant 
one). Even if the fli are identified, finding all gi is gener- 
ally impossible, unless we know that one dominates, and 
we know which one that is. The only simple case is if all 
phonons have roughly equal frequencies, when one can 
treat them as a single mode with coupling g^ s = 9i • 

To summarize, we have found a generalization of the 
simple, yet accurate Momentum Average approximations 
to the problem of Holstein-type coupling to multiple 
phonon modes. Our results show that even perturba- 
tionally weak coupling to a second phonon can lead to 
essential changes of the spectral weight, if its frequency 
is less than that of the dominant phonon. In such cases, 
the simple way of extracting the electron-phonon cou- 
pling from ARPES data is likely to lead to wrong values. 
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